home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dlaexc.f < prev    next >
Text File  |  1996-07-19  |  11KB  |  356 lines
  1.       SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
  2.      $                   INFO )
  3. *
  4. *  -- LAPACK auxiliary routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     February 29, 1992
  8. *
  9. *     .. Scalar Arguments ..
  10.       LOGICAL            WANTQ
  11.       INTEGER            INFO, J1, LDQ, LDT, N, N1, N2
  12. *     ..
  13. *     .. Array Arguments ..
  14.       DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WORK( * )
  15. *     ..
  16. *
  17. *  Purpose
  18. *  =======
  19. *
  20. *  DLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
  21. *  an upper quasi-triangular matrix T by an orthogonal similarity
  22. *  transformation.
  23. *
  24. *  T must be in Schur canonical form, that is, block upper triangular
  25. *  with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block
  26. *  has its diagonal elemnts equal and its off-diagonal elements of
  27. *  opposite sign.
  28. *
  29. *  Arguments
  30. *  =========
  31. *
  32. *  WANTQ   (input) LOGICAL
  33. *          = .TRUE. : accumulate the transformation in the matrix Q;
  34. *          = .FALSE.: do not accumulate the transformation.
  35. *
  36. *  N       (input) INTEGER
  37. *          The order of the matrix T. N >= 0.
  38. *
  39. *  T       (input/output) DOUBLE PRECISION array, dimension (LDT,N)
  40. *          On entry, the upper quasi-triangular matrix T, in Schur
  41. *          canonical form.
  42. *          On exit, the updated matrix T, again in Schur canonical form.
  43. *
  44. *  LDT     (input)  INTEGER
  45. *          The leading dimension of the array T. LDT >= max(1,N).
  46. *
  47. *  Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
  48. *          On entry, if WANTQ is .TRUE., the orthogonal matrix Q.
  49. *          On exit, if WANTQ is .TRUE., the updated matrix Q.
  50. *          If WANTQ is .FALSE., Q is not referenced.
  51. *
  52. *  LDQ     (input) INTEGER
  53. *          The leading dimension of the array Q.
  54. *          LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.
  55. *
  56. *  J1      (input) INTEGER
  57. *          The index of the first row of the first block T11.
  58. *
  59. *  N1      (input) INTEGER
  60. *          The order of the first block T11. N1 = 0, 1 or 2.
  61. *
  62. *  N2      (input) INTEGER
  63. *          The order of the second block T22. N2 = 0, 1 or 2.
  64. *
  65. *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
  66. *
  67. *  INFO    (output) INTEGER
  68. *          = 0: successful exit
  69. *          = 1: the transformed matrix T would be too far from Schur
  70. *               form; the blocks are not swapped and T and Q are
  71. *               unchanged.
  72. *
  73. *  =====================================================================
  74. *
  75. *     .. Parameters ..
  76.       DOUBLE PRECISION   ZERO, ONE
  77.       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  78.       DOUBLE PRECISION   TEN
  79.       PARAMETER          ( TEN = 1.0D+1 )
  80.       INTEGER            LDD, LDX
  81.       PARAMETER          ( LDD = 4, LDX = 2 )
  82. *     ..
  83. *     .. Local Scalars ..
  84.       INTEGER            IERR, J2, J3, J4, K, ND
  85.       DOUBLE PRECISION   CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
  86.      $                   T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2,
  87.      $                   WR1, WR2, XNORM
  88. *     ..
  89. *     .. Local Arrays ..
  90.       DOUBLE PRECISION   D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
  91.      $                   X( LDX, 2 )
  92. *     ..
  93. *     .. External Functions ..
  94.       DOUBLE PRECISION   DLAMCH, DLANGE
  95.       EXTERNAL           DLAMCH, DLANGE
  96. *     ..
  97. *     .. External Subroutines ..
  98.       EXTERNAL           DLACPY, DLANV2, DLARFG, DLARFX, DLARTG, DLASY2,
  99.      $                   DROT
  100. *     ..
  101. *     .. Intrinsic Functions ..
  102.       INTRINSIC          ABS, MAX
  103. *     ..
  104. *     .. Executable Statements ..
  105. *
  106.       INFO = 0
  107. *
  108. *     Quick return if possible
  109. *
  110.       IF( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 )
  111.      $   RETURN
  112.       IF( J1+N1.GT.N )
  113.      $   RETURN
  114. *
  115.       J2 = J1 + 1
  116.       J3 = J1 + 2
  117.       J4 = J1 + 3
  118. *
  119.       IF( N1.EQ.1 .AND. N2.EQ.1 ) THEN
  120. *
  121. *        Swap two 1-by-1 blocks.
  122. *
  123.          T11 = T( J1, J1 )
  124.          T22 = T( J2, J2 )
  125. *
  126. *        Determine the transformation to perform the interchange.
  127. *
  128.          CALL DLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP )
  129. *
  130. *        Apply transformation to the matrix T.
  131. *
  132.          IF( J3.LE.N )
  133.      $      CALL DROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS,
  134.      $                 SN )
  135.          CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
  136. *
  137.          T( J1, J1 ) = T22
  138.          T( J2, J2 ) = T11
  139. *
  140.          IF( WANTQ ) THEN
  141. *
  142. *           Accumulate transformation in the matrix Q.
  143. *
  144.             CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
  145.          END IF
  146. *
  147.       ELSE
  148. *
  149. *        Swapping involves at least one 2-by-2 block.
  150. *
  151. *        Copy the diagonal block of order N1+N2 to the local array D
  152. *        and compute its norm.
  153. *
  154.          ND = N1 + N2
  155.          CALL DLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD )
  156.          DNORM = DLANGE( 'Max', ND, ND, D, LDD, WORK )
  157. *
  158. *        Compute machine-dependent threshold for test for accepting
  159. *        swap.
  160. *
  161.          EPS = DLAMCH( 'P' )
  162.          SMLNUM = DLAMCH( 'S' ) / EPS
  163.          THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
  164. *
  165. *        Solve T11*X - X*T22 = scale*T12 for X.
  166. *
  167.          CALL DLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD,
  168.      $                D( N1+1, N1+1 ), LDD, D( 1, N1+1 ), LDD, SCALE, X,
  169.      $                LDX, XNORM, IERR )
  170. *
  171. *        Swap the adjacent diagonal blocks.
  172. *
  173.          K = N1 + N1 + N2 - 3
  174.          GO TO ( 10, 20, 30 )K
  175. *
  176.    10    CONTINUE
  177. *
  178. *        N1 = 1, N2 = 2: generate elementary reflector H so that:
  179. *
  180. *        ( scale, X11, X12 ) H = ( 0, 0, * )
  181. *
  182.          U( 1 ) = SCALE
  183.          U( 2 ) = X( 1, 1 )
  184.          U( 3 ) = X( 1, 2 )
  185.          CALL DLARFG( 3, U( 3 ), U, 1, TAU )
  186.          U( 3 ) = ONE
  187.          T11 = T( J1, J1 )
  188. *
  189. *        Perform swap provisionally on diagonal block in D.
  190. *
  191.          CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
  192.          CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
  193. *
  194. *        Test whether to reject swap.
  195. *
  196.          IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 3,
  197.      $       3 )-T11 ) ).GT.THRESH )GO TO 50
  198. *
  199. *        Accept swap: apply transformation to the entire matrix T.
  200. *
  201.          CALL DLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK )
  202.          CALL DLARFX( 'R', J2, 3, U, TAU, T( 1, J1 ), LDT, WORK )
  203. *
  204.          T( J3, J1 ) = ZERO
  205.          T( J3, J2 ) = ZERO
  206.          T( J3, J3 ) = T11
  207. *
  208.          IF( WANTQ ) THEN
  209. *
  210. *           Accumulate transformation in the matrix Q.
  211. *
  212.             CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
  213.          END IF
  214.          GO TO 40
  215. *
  216.    20    CONTINUE
  217. *
  218. *        N1 = 2, N2 = 1: generate elementary reflector H so that:
  219. *
  220. *        H (  -X11 ) = ( * )
  221. *          (  -X21 ) = ( 0 )
  222. *          ( scale ) = ( 0 )
  223. *
  224.          U( 1 ) = -X( 1, 1 )
  225.          U( 2 ) = -X( 2, 1 )
  226.          U( 3 ) = SCALE
  227.          CALL DLARFG( 3, U( 1 ), U( 2 ), 1, TAU )
  228.          U( 1 ) = ONE
  229.          T33 = T( J3, J3 )
  230. *
  231. *        Perform swap provisionally on diagonal block in D.
  232. *
  233.          CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
  234.          CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
  235. *
  236. *        Test whether to reject swap.
  237. *
  238.          IF( MAX( ABS( D( 2, 1 ) ), ABS( D( 3, 1 ) ), ABS( D( 1,
  239.      $       1 )-T33 ) ).GT.THRESH )GO TO 50
  240. *
  241. *        Accept swap: apply transformation to the entire matrix T.
  242. *
  243.          CALL DLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK )
  244.          CALL DLARFX( 'L', 3, N-J1, U, TAU, T( J1, J2 ), LDT, WORK )
  245. *
  246.          T( J1, J1 ) = T33
  247.          T( J2, J1 ) = ZERO
  248.          T( J3, J1 ) = ZERO
  249. *
  250.          IF( WANTQ ) THEN
  251. *
  252. *           Accumulate transformation in the matrix Q.
  253. *
  254.             CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
  255.          END IF
  256.          GO TO 40
  257. *
  258.    30    CONTINUE
  259. *
  260. *        N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
  261. *        that:
  262. *
  263. *        H(2) H(1) (  -X11  -X12 ) = (  *  * )
  264. *                  (  -X21  -X22 )   (  0  * )
  265. *                  ( scale    0  )   (  0  0 )
  266. *                  (    0  scale )   (  0  0 )
  267. *
  268.          U1( 1 ) = -X( 1, 1 )
  269.          U1( 2 ) = -X( 2, 1 )
  270.          U1( 3 ) = SCALE
  271.          CALL DLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 )
  272.          U1( 1 ) = ONE
  273. *
  274.          TEMP = -TAU1*( X( 1, 2 )+U1( 2 )*X( 2, 2 ) )
  275.          U2( 1 ) = -TEMP*U1( 2 ) - X( 2, 2 )
  276.          U2( 2 ) = -TEMP*U1( 3 )
  277.          U2( 3 ) = SCALE
  278.          CALL DLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 )
  279.          U2( 1 ) = ONE
  280. *
  281. *        Perform swap provisionally on diagonal block in D.
  282. *
  283.          CALL DLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK )
  284.          CALL DLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK )
  285.          CALL DLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK )
  286.          CALL DLARFX( 'R', 4, 3, U2, TAU2, D( 1, 2 ), LDD, WORK )
  287. *
  288. *        Test whether to reject swap.
  289. *
  290.          IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ),
  291.      $       ABS( D( 4, 2 ) ) ).GT.THRESH )GO TO 50
  292. *
  293. *        Accept swap: apply transformation to the entire matrix T.
  294. *
  295.          CALL DLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK )
  296.          CALL DLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK )
  297.          CALL DLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK )
  298.          CALL DLARFX( 'R', J4, 3, U2, TAU2, T( 1, J2 ), LDT, WORK )
  299. *
  300.          T( J3, J1 ) = ZERO
  301.          T( J3, J2 ) = ZERO
  302.          T( J4, J1 ) = ZERO
  303.          T( J4, J2 ) = ZERO
  304. *
  305.          IF( WANTQ ) THEN
  306. *
  307. *           Accumulate transformation in the matrix Q.
  308. *
  309.             CALL DLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK )
  310.             CALL DLARFX( 'R', N, 3, U2, TAU2, Q( 1, J2 ), LDQ, WORK )
  311.          END IF
  312. *
  313.    40    CONTINUE
  314. *
  315.          IF( N2.EQ.2 ) THEN
  316. *
  317. *           Standardize new 2-by-2 block T11
  318. *
  319.             CALL DLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ),
  320.      $                   T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN )
  321.             CALL DROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT,
  322.      $                 CS, SN )
  323.             CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
  324.             IF( WANTQ )
  325.      $         CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
  326.          END IF
  327. *
  328.          IF( N1.EQ.2 ) THEN
  329. *
  330. *           Standardize new 2-by-2 block T22
  331. *
  332.             J3 = J1 + N2
  333.             J4 = J3 + 1
  334.             CALL DLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ),
  335.      $                   T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN )
  336.             IF( J3+2.LE.N )
  337.      $         CALL DROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ),
  338.      $                    LDT, CS, SN )
  339.             CALL DROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN )
  340.             IF( WANTQ )
  341.      $         CALL DROT( N, Q( 1, J3 ), 1, Q( 1, J4 ), 1, CS, SN )
  342.          END IF
  343. *
  344.       END IF
  345.       RETURN
  346. *
  347. *     Exit with INFO = 1 if swap was rejected.
  348. *
  349.    50 CONTINUE
  350.       INFO = 1
  351.       RETURN
  352. *
  353. *     End of DLAEXC
  354. *
  355.       END
  356.